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ABSTRACT* 

It  is  shown  that  for  quasi-1 inear  hyperbolic  systems  of  the  conservation 
form  *  -  Fx  =  -  AWx>  it  is  possible  to  build  up  relatively  simple  finite 
difference  numerical  schemes  accurate  to  3rd  and  4th  order  provided  that  the 
matrix  A  satisfies  commutativity  relations  with  its  partial-derivative-matrices. 
This  requirement  is  not  fulfilled  by  any  known  physical  systems  of  equations. 

These  schemes  generalize  the  Lax-Wendroff  2nd  order  one,  and  are  written  down 
explicitly.  As  found  by  Strang,  odd  order  schemes  are  linearly  unstable  unless 
modified  by  adding  a  term  containing  the  next  higher  space  derivative.  Thus 
stabilized,  the  schemes,  both  odd  and  even,  can  be  made  to  meet  the  C.F.L 
(Courant-Friedrichs-Lewy)  criterion.  Numerical  calculations  were  made  with  a 
3rd  order  and  a  4th  order  scheme  for  scalar  equations  with  continuous  and  dis¬ 
continuous  solutions.  The  results  are  compared  with  analytic  solutions  and  the 
predicted  improvement  is  verified. 

The  computation  reported  herein,  were  carried  out  on  the  CDC-3400 
computer  at  the  Tel  Aviv  University  computation  center. 
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1.  INTRODUCTION 


When  dealing  with  one-dimensional  problems  in  continuum  mechanics,  and 
in  particular  hydrodynamics,  it  is  often  necessary  to  solve  non-linear  hyper¬ 
bolic  systems  of  the  form 


where  (  and  (  )  denote,  respectively,  partial  differentiation  wi th 
respect  to  the  time  and  space  coordinates.  W  is  a  vector  whose  components 
are  the  unknown  functions  and  F  is  a  vector  whose  components  are  dependent 
functionally  on  the  components  of  W  only.  We  consider  the  quasi-linear  case 
wherein  F^  =  AW^,  A  being  a  matrix  whose  components  depend  on  the  unknown 
functions  only,  and  not  on  their  derivatives.  Since  eq.  (1)  is  assumed  to  be 
hyperbolic,  the  eigenvalues  of  A  are  all  real. 

Systems  of  the  form  of  equation  (l)  are  called  "Conservation  Law  Form" 
systems.  Various  numerical  schemes  for  their  solutions  have  been  developed 
[1],  [2],  [3],  [4],  [ 5 1  starting  with  Lax  and  Wendroff  [1], 


Keeping  in  mind  that  ultimately  the  main  interest  will  focus  on  multi¬ 
dimensional  systems  (say  x,y,z  and  t  or  x,r,$  and  t)  it  is  of  obvious 
importance  to  develop  numerical  schemes  whose  order  of  accuracy  is  higher 
than  one.  A  widely  used  2nd  order  accuracy  scheme  is  the  one  due  to  Lax  and 
Wendroff  [1].  Their  finite  differences  approximation  is  written  thusly: 
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(n  =  0, 1 ,2, _  ) 

!f  the  problem  contains  discontinuities  (such  as  shocks  which  may  develop 
even  if  the  initial  conditions  are  smooth;  see  [6])  the  system  may  be  handled 
either  by  adding  a  non-linear  artificial  viscosity  term  [1]  or  by  iterative 
methods  [4],  A  stability  criterion  determines  Atn  for  the  predetermined  and 
fixed  Ax. 

With  view  in  mind  towards  multi-dimensional  computations,  it  is  of 
interest  to  consider  3rd  order  accuracy  schemes  for  non-linear  hyperbolic 
equations  of  the  type  (1)  where  for  the  present  we  shall  consider  the  scalar 
case,  i.e.  a  simple  conservation-form  equation.  A  first  attempt  in  this  direction 
is  due  to  Burstein  [7]  who  developed  a  three-step  approach  in  analogy  to 
P.ichtmyer's  two-step  method  [3]  which  approximates  the  Lax-Wendroff  2nd  order 
scheme.  We  shall  use  the  basic  ideas  of  Lax  and  Wendroff  [1]  for  estimating 
truncation  errors  in  order  to  construct,  a  third  order  scheme. 


2.  DERIVATION  OF  THE  METHOD 


The  Lax-Wendroff  (LW)  method  is  based  on  the  fact  that  from  the  equation 


1 
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This  allows  the  construction  of  a  2nd  order  scheme  by  developing  W(x,t,  +  At) 
in  a  Taylor  series  which  to  order  (At)2  is  given  by 


W(x,t  +  At)  =  W(x,t)  +  (At)Wt  +  +  0(At3) 


and  the  time  derivatives  are  replaced  by  space  derivatives  through  the  use  of 
equations  (1)  and  (3).  This  allows  the  build  up  of  a  finite  differences  scheme 
which  steps  up  W  in  time  by  using  only  spatial  differences. 

Our  first  task  is  to  construct  a  formula  similar  to  equation  (3)  for 
higher  order  time  derivatives.  We 'make  the  following  claim: 

If  the  matrix  A  of  the  hyperbolic  system  (1)  is  commutative 
with  its  partial-derivative-matrices  then 


s!!“ .  (-0"  iLL(fln-,f ) 
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for  every  natural  number  n. 


The  proof  proceeds  as  follows: 


W  =  (AF  )  =  A  F  +  AF 

tt  xx  xx  xx 


on  the  other  hand. 
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Comparing  (5)  and  (6): 


AW  =  -A  F  =  -A  AW 
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With  these  preliminaries  and  with  our  assertion  (4)  known  to  be  true  for 
n  =  1  and  2  (even  without  the  commutativity  condition),  we  now  consider 
the  case  n  =  3: 


“ttt  -  -  I(AFx>tix  -  [\Fx  *  AFxt’x  ‘  lAtA“x  +  A<-AFx>x1x-  <8) 

Now  substitute  the  commutativity  restriction  A^A  =  AA{  into  (8)  to  get 

W  «  [AA.W  -  A(A  F  +  AF  )] 
ttt  t  x  xx  xx  x 

Using  (7)  we  obtain  (using  A^A  =  AA^): 
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We  have  thus  verified  our  claim  (eq.{4))  for  the  case  n  =  3-  i*  ;s  easy  to 
show  by  induction  that  equation  (4)  is  valid  for  all  n. 


With  the  aid  of  this  result,  we  can  construct  a  finite  difference  scheme 
to  any  desired  truncation-error  by  writing  the  required  Taylor  series: 

W(x, t  +  At)  =  W(x,0  -  7  (")k  ^l-4rr(Ak',F  )  +  OUAt)1^1]  .  (10) 

k=l  k!  3x 

One  has  only  to  take  care  to  represent  the  various  derivatives  by  a  finite 
differences  expression  which  has  the  proper  accuracy  so  that  the  overall 
scheme  retains  the  desired  truncation  error.  It  should  be  noted  that  by  the 
Cayley-Hami 1  ton  Tht orem  it  is  possible  to  express  A  (k  i  r)  in  terms  of 
A,  A2,...,  Ar  *  where  A  is  of  order  r  x  r. 


UJ 


A  THiRD  ORDER  FINITE  DIFFERENCES  SCHEME 


We  shall  consider  now  the  case  of  a  system  obeying  the  restriction 
on  A: 

W  «  F  =  A(W)W  (11) 

t  x  x 

with  the  initial  condition 

W(x,o)  =  $(x)  .  (12) 

In  order  for  all  the  terms  in  the  finite  differences  representation  to  be  of 
at  least  of  third  order,  it  is  necessary  to  improve  the  representation  of 
the  first  derivative  thusly: 


It  should  be  noted  that  for  a  third  order  accuracy  scheme  one  may  still  use 


fu"  ♦  w"l 

Aj+|/2  *  AM - ^ or  something  else  equivalent;  but  for  hi 

schemes  it  will  be  necessary  to  utilize  a  better  interpolation. 


or  something  else  equivalent;  but  for  higher  order 


We  can  now  write  down  a  finite  differences  scheme  of  third  order 
accuracy: 


“T1  *  *  «Ki  -  7>>  -  -  >7,  *  *7,  -  f;.2» 
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with  6  -  e  -  I.  If  we  set  6  =  s  =  0  we  get  back  to  the  2nd  order  accuracy 
scheme  of  la*  and  Wendroff.  As  shall  see,  the  scheme  (Ik)  Is  unstable  as 
It  stands,  and  in  order  to  stabilize  It  we  have  to  add  artificial  viscosity 
terms.  This  state  of  affairs  is  typical  in  schemes  of  odd  order  of  accuracy. 


4.  THE  STABILIZING  ARTIFICIAL  VISCOSITY 


We  propose  an  artificial  viscosity  term  of  the  form 


2jj-h4(AzA2  +  vA4A4)W^ 


.  34W* 


*4*  - J 

**  3x4 


The  part  of  (15)  which  is  proportional  to  \*A\x  is  due  to  the  expressing 

of  to  a  higher  accuracy  than  0(h2) .  This  improvement,  i f  given  in  a 

conservation  form,  is  written  thusly: 
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It  is  the  second  term  on  the  R.H.S.  of  ( 1 6)  which  gives  rise  to  \2A2W^x 
(with  A  being  taken  constant,  since  the  term  is  not  needed  for  accuracy, 
only  for  stabi I i ty) . 


The  a4A4W^  part  of  (15),  however,  is  due  to  the  extending  of  the 
Taylor  series;  i.e.  it  comes  from 


V.,  =  (A3F  ) 

At  x  xxx 


(17) 


The  finite  difference  representation  of  (17)  is,  to  3rd  order  accuracy. 
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In  practice  we'll  take  the  artificial  viscosity  term  either  in  its  conservation 
form  (i.e-  use  ( 1 6)  and  (18),  or  in  the  linear  version  where  the  finite 
difference  form  of  (15)  is: 


#>-2(a^): 


+  vA4(a")4][W^2  -  Wj+J  +  6w”  -  4W"_,  +  w"_2] 


(19) 


If  we  examine  the  linear  stability  of  the  scheme  in  the  usual  Von-"eumann 
fashion,  then  we  find  that  fo-  a*e»6*-v«l  the  criterion  is  Aa  $  1 
where  a  is  the  spectral  radius  of  A.  This  result  is  also  a  special  case  of 
Theorem  1  in  Strang's  paper  of  1962.(8], 

In  principal,  it  is  possible  to  build  up  in  the  above  manner,  numerical 
schemes  of  any  desired  accuracy.  It  seems  that  such  schemes  of  odd  orders  will 
not  be  linearly  stable  unless  an  artificial  viscosity  term  is  added.  The 
artificial  viscosity  term  that  we  have  added  contains  a  term  which  is 
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proportional  to  the  next  higher  even  derivative  of  the  Taylor  series 
development  of  W(x,t  +  At).  Thus  it  might  be  possible  to  take  just  that 
amount  of  artificial  viscosity  which  will  not  only  stabilize  the  scheme,  but 
also  will  add  one  more  order  of  truncation  accuracy.  Of  course,  when  doing 
this,  care  has  to  be  taken  that  all  the  differencing  is  consistent  with  the 
higher  order  accuracy.  In  effect,  this  is  the  case  with  the  Lax-Wendroff 


term,  which  might  be  considered  as  a  stabilizing  term  added  to  a  first  order 

A 

scheme.  If  the  coefficient  of  (AF  )  is  cleverly  chosen  to  be  yp*.  then 
an  additional  bonus  is  the  2nd  order  accuracy. 
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We'll  take  these  solutions  at  x  =  0  and  x  =  1  to  serve  as  the  boundary 
conditions  for  the  numerical  work  which  is  to  be  checked  in  0  <  x  <  1 
against  the  above  analytic  solutions. 

The  case  of  a  solution  containing  a  discontinuity  which  is  created  at 
some  t  =  t  is  demonstrated  by  taking  the  following  initial  conditions: 


»(x) 


-  »  <  X  £  0 


$  x  s  2© 


20  $  x  < 


The  solution  is  given  by: 


u(x,t)  = 


2e  -  x 
0  -  t 


-  »  <  x  £  t+0 


t  +  0  <  x  r  25 


20  $  x  <  » 


(0  <  t  <  S) 


and  for  t  2  t  =  O 


-  «  <  x  f 


t  +  30 


u(x,t)  = 


(t  *  0) 


t  +  33 


<  X  $  1 


The  "shock-wave”  is  created  at  time  t  =  0  end  at  the  location  x  —  20, 
and  moves  to  the  right  with  the  speed  xg  =  1/2.  In  the  numerical  com¬ 
putations  corresponding  to  this  case,  we  shall  examine  the  behavior  of  the 
solution  also  across  the  discontinuity  and  compare  it  to  both  the  analytic 
solution  and  the  results  obtained  from  the  standard  Lax-Wendroff  method. 


6,  NUMERICAL  RESULTS 


In  reporting  on  the  numerical  work,  we  shall  compare  the  standard 
Lax-Wendroff  results  with  those  of  our  third  order  scheme  (eq.  (1*0), 
either  with  the  linear  viscosity  (eq.  (19))  or  with  the  conservation  form 
artificial  viscosity  ( ( 1 6)  +  (18)). 


6a.  AN  EXAMPLE  FOR  SMOOTH  SOLUTIONS 


We  took  Ax  =  0.005,  u(x,o)  =  $(x)  =  x2  and,  from  (22), 
u(o,t)  =  0 

/,  .i  2t  +  l  -  /iTi 

u  ( 1 ,  t)  «  - 

2t2 

As  expected,  the  maximum  relative  error  using  the  scheme  presented  above 
is  about  100  times  smaller  than  the  one  given  by  the  standard  Lax-Wendrcff 
method.  This  is  the  case  when  we  use  the  linear  artificial  viscosity 
(eq.  (19)).  When  the  conservation  form  of  the  artificial  viscosity  is  used, 
the  ratio  of  the  maximum  relative  errors  decreases  from  about  1/100  to 
about  1/1000. 
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6b.  AN  EXAMPLE  FOR  DISCONTINUOUS  SOLUTIONS 


Here  we  took  +  uu^  =0  with  the  initial  distribution  (25)  and 
boundary  conditions  according  to  (26)  and  (27).  The  results  indicate  that 
the  3rd  order  accuracy  scheme  gives  a  slight  improvement  over  the  L-W 
calculation  in  the  large  gradient  region,  in  the  sense  that  the  ‘shock* 
is  slightly  steeper  and  the  post-shock  oscillations  are  weaker  and  are 
damped  more  quickly.  On  the  other  hand,  unlike  the  L-W  case,  there  is  a 
very  small  negative  perturbation  ahead  of  the  "shock".  The  appearance  of 
this  precursor  perturbation  is  due  to  the  fact  that  in  the  3rd  order 
scheme  one  uses,  in  addition  to  u^  and  uj+j»  also  uj+2* 


7.  A  FOURTH  ORDER  FINITE  DIFFERENCES  SCHEME 


As  was  mentioned  before,  if  in  the  term  stemming  from  we  represent 

^*■+1/2  a  i gfier  order  interpolation  formula,  and  if  we  take  *2=1,  then 
our  scheme  becomes  of  fourth  order  accuracy.  This  is  in  line  with  the  remarks 
at  the  end  of  Section  3  -  adding  a  stabilizing  term  to  an  odd  order  scheme 
can  raise  the  order  cf  accuracy  ?f  the  coefficients  of  this  added  derivative 
are  chosen  properly. 


In  the  present  case,  the  fourth  order  finite  difference  scheme  has 
the  following  form: 
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where 


A°  =  — (An  +  A°)  -  — »-{An  -  A°  -  a°  +  A°  ) 

j+l/2  rjil  j  T5j±2  j±1  j  j-r 


The  scheme  (28)  meets  the  C.F.L  stability  criterion.  Note  that  if  i»  *  0 
then  eq.  (28)  is  the  3rd  order  scheme  with  the  linear  artificial  viscosity 
represented  in  a  conservative  form.  If  also  5  =  e  =  a  =  0,  then  we  have 
the  Lax-.*«ndroff  scheme.  For  the  fourth  order  accuracy,  we  must  use 
5  =  e  =  Q  =  ii  =  l.  We  ran  test  runs  with  $(x)  =  xa  in  order  to  determine 
in  practice,  and  compare  to  prediction,  the  amount  by  which  the  grid  can  be 
coarsened  and  still  maintain  the  same  maximum  relative  error  as  we  go  from 
Lax-Wendroff  to  3rd  order  and  then  4th  order  schemes.  The  grid  sizes  (based 


iAiiUtVvUi 
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on  Ax  =•  .005  for  the  L-W  scheme)  were  found  to  be,  respectively: 

Ax,  ,,  =  .1/200;  Ax,  ,  ,  =  .1/50;  and  Ax,...  .  =  .1/25.  These 

L.W.  3rd  order  ’  _6  *»th  order 

grids  produce  absolute  errors  of  C*10  where  for  o  =  2  and  t  =  3.» 

1  <  C  <  5. 

In  conclusion,  it  may  be  stated  that  a  3rd  or  ^th  order  scheme,  such 
as  the  ones  proposed  in  this  paper,  will  in  the  smooth  part  of  the  solution 
to  a  hyperbolic  problem,  yield  in  practice  as  well  as  in  theory,  one  or  two 
orders  of  accuracy  higher  than,  say,  the  Lax-Wendroff  method.  In  the  regions 
near  shock-like  discontinuities,  the  improvement  over  the  L-W  scheme  is  not 
as  great.  This  opens  up  the  possibility  that  in  multidimensional  cases  we 
shall  be  able  in  a  practical  manner  to  overcome  partially  the  problems  of 
restricted  machine  memory  by  using  coarser  grids  and  higher  order  schemes. 
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